#include "wrapper_hypre.h"

typedef HYPRE_Int idx_t;
typedef HYPRE_Real data_t;

const idx_t num_diag = 7;

int main(int argc, char* argv[])
{
    int my_pid, num_procs;
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &my_pid);
    MPI_Comm_size(MPI_COMM_WORLD, &num_procs);

    idx_t arg_cnt = 1;
    const idx_t local_m = atoi(argv[arg_cnt++]),// 每个part内有一个box，该box有m*m*m个cell
                np_per_dim = atoi(argv[arg_cnt++]),// 每个part内有一个box，该box有p*p*p个进程组成
                nparts = 2;
    const std::string solver_name = std::string(argv[arg_cnt++]),
                        prec_name = std::string(argv[arg_cnt++]);
    const idx_t ncells_per_dim = local_m * np_per_dim,
                ncells_per_part = ncells_per_dim * ncells_per_dim * ncells_per_dim,
                nprocs_per_part = np_per_dim * np_per_dim * np_per_dim;
    assert(ncells_per_dim % 4 == 0);
    assert(num_procs == nprocs_per_part * nparts);

    const idx_t my_part = my_pid / nprocs_per_part,
                pid_in_part = my_pid - my_part * nprocs_per_part;
    idx_t cart_ids[3];
    {// 计算本进程在part内处于什么位置
        cart_ids[0] = pid_in_part / (np_per_dim * np_per_dim);
        cart_ids[1] = (pid_in_part - cart_ids[0] * np_per_dim * np_per_dim) / np_per_dim;
        cart_ids[2] = (pid_in_part - cart_ids[0] * np_per_dim * np_per_dim - cart_ids[1] * np_per_dim);
    }

    idx_t my_ilower[3], my_iupper[3];
    for (idx_t d = 0; d < 3; d++) {
        my_ilower[d] =  cart_ids[d]     * local_m;
        my_iupper[d] = (cart_ids[d] + 1)* local_m - 1;
    }
#ifdef _zyDEBUG
    for (int p = 0; p < num_procs; p++) {
        if (my_pid == p) {
            printf("Proc %d part %d: (%d,%d,%d)~(%d,%d,%d)\n", my_pid, my_part,
                my_ilower[0], my_ilower[1], my_ilower[2],
                my_iupper[0], my_iupper[1], my_iupper[2] );
            fflush(stdout);
        }
        MPI_Barrier(MPI_COMM_WORLD);
    }
    MPI_Barrier(MPI_COMM_WORLD);
#endif
    idx_t obj_type = HYPRE_SSTRUCT;
    if (prec_name == "AMG") obj_type = HYPRE_PARCSR;
    if (my_pid == 0) {
        printf("\033[1;35m Local M = %d P: %d Total #procs: %d Total #dof: %d \033[0m \n", local_m, np_per_dim, num_procs, ncells_per_part * nparts);
        printf("object type : %d\n", obj_type);
    }

    HYPRE_Init();

    HYPRE_SStructGrid ssgrid;
    HYPRE_SStructStencil stencil;
    HYPRE_SStructGraph ssgraph;
    HYPRE_SStructMatrix   A;
    HYPRE_SStructVector   b, x, y, x0;
    HYPRE_ParCSRMatrix par_A;
    HYPRE_ParVector par_b, par_x, par_y;
    TEST_RECORD records[TEST_CNT];
    TEST_CONFIG config;
    config.case_name = "SAMR_bench";
    config.rtol = 1.0e-7;

    {
        HYPRE_SStructGridCreate(MPI_COMM_WORLD, 3, nparts, &ssgrid);
        HYPRE_SStructGridSetExtents(ssgrid, my_part, my_ilower, my_iupper);
        HYPRE_SStructVariable vtypes[1] = { HYPRE_SSTRUCT_VARIABLE_CELL };
        for (idx_t part = 0; part < nparts; part++) {// 每个part都要执行定义变量的操作
            HYPRE_SStructGridSetVariables(ssgrid, part, sizeof(vtypes)/sizeof(HYPRE_SStructVariable), vtypes);
        }

        HYPRE_SStructGridAssemble(ssgrid);// a collective call finalizing the grid assembly.

        HYPRE_SStructStencilCreate(3, num_diag, &stencil);
        idx_t offsets[7][3] = {
            {-1,0,0}, {0,-1,0}, {0,0,-1}, {0,0,0}, {0,0,1}, {0,1,0}, {1,0,0},
        };

        for (idx_t e = 0; e < num_diag; e++) {
            HYPRE_SStructStencilSetEntry(stencil, e, offsets[e], 0);// 只有0号变量一种
        }
        
        HYPRE_SStructGraphCreate(MPI_COMM_WORLD, ssgrid, &ssgraph);
        HYPRE_SStructGraphSetObjectType(ssgraph, obj_type);
        for (idx_t part = 0; part < nparts; part++) {// 每个part都要执行定义模板的操作
            HYPRE_SStructGraphSetStencil(ssgraph, part, 0, stencil);
        }

        // 添加上非stencil的连接关系
        const idx_t amr_stride = 2;
        assert(amr_stride == 2);
        const idx_t c_lbdr = ncells_per_dim / 4     - 1,// 原网格中边界处较下者
                    c_ubdr = ncells_per_dim / 4 * 3    ;// 原网格中边界处较上者
        if (my_pid == 0) printf("lbdr/ubdr @ C: %d %d\n", c_lbdr, c_ubdr);
        if (my_part == 0) {// 原来的网格
            idx_t ngb_part = 1;
            for (idx_t d = 2; d >= 0; d--) {// 在三个方向逐次添加连接关系
                if (my_ilower[d] <= c_lbdr && c_lbdr <= my_iupper[d]) {
                    idx_t bdr_ilower[3], bdr_iupper[3];// 确定原网格上的边界位置
                    for (idx_t _k = 0; _k < 3; _k++) {
                        if (_k == d) {
                            bdr_ilower[_k] = c_lbdr;
                            bdr_iupper[_k] = c_lbdr;
                        } else {
                            bdr_ilower[_k] = std::max(c_lbdr + 1, my_ilower[_k]);
                            bdr_iupper[_k] = std::min(c_ubdr - 1, my_iupper[_k]);
                        }
                    }
#ifdef _zyDEBUG
                    printf("Proc %d bdr (%d,%d,%d) ~ (%d,%d,%d)\n", my_pid,
                        bdr_ilower[0], bdr_ilower[1], bdr_ilower[2],
                        bdr_iupper[0], bdr_iupper[1], bdr_iupper[2] );
#endif
                    // 从粗往细看，自己的一个粗点对应4个细点
                    for (idx_t fngb = 0; fngb < 4; fngb ++) {
                        idx_t offset[3];
                        if (d == 0) {
                            offset[0] = __builtin_popcount(0          ); offset[1] = __builtin_popcount(fngb & 0b01); offset[2] = __builtin_popcount(fngb & 0b10);
                        } else if (d == 1) {
                            offset[0] = __builtin_popcount(fngb & 0b01); offset[1] = __builtin_popcount(0          ); offset[2] = __builtin_popcount(fngb & 0b10);
                        } else if (d == 2) {
                            offset[0] = __builtin_popcount(fngb & 0b01); offset[1] = __builtin_popcount(fngb & 0b10); offset[2] = __builtin_popcount(0          );
                        }
                        for (idx_t i2 = bdr_ilower[2]; i2 <= bdr_iupper[2]; i2++)
                        for (idx_t i1 = bdr_ilower[1]; i1 <= bdr_iupper[1]; i1++)
                        for (idx_t i0 = bdr_ilower[0]; i0 <= bdr_iupper[0]; i0++) {
                            idx_t index[3] = {i0, i1, i2}, to_index[3];
                            for (idx_t _k = 0; _k < 3; _k++) {
                                if (_k == d) to_index[_k] = 0;// 注意此时对应加密网格的下边界
                                else         to_index[_k] = (index[_k] - (c_lbdr + 1)) * amr_stride + offset[_k];
                            }
#ifdef _zyDEBUG
                            hypre_printf("index: [%d](%d, %d, %d) - to_index: [%d](%d, %d, %d)\n",
                                        my_part, index[0], index[1], index[2],
                                        ngb_part, to_index[0], to_index[1], to_index[2]);
#endif
                            HYPRE_SStructGraphAddEntries(ssgraph, my_part, index, 0, ngb_part, to_index, 0);
                        }
                    }
                }// if lower bdr included
                if (my_ilower[d] <= c_ubdr && c_ubdr <= my_iupper[d]) {
                    idx_t bdr_ilower[3], bdr_iupper[3];
                    for (idx_t _k = 0; _k < 3; _k++) {
                        if (_k == d) {
                            bdr_ilower[_k] = c_ubdr;
                            bdr_iupper[_k] = c_ubdr;
                        } else {
                            bdr_ilower[_k] = std::max(c_lbdr + 1, my_ilower[_k]);
                            bdr_iupper[_k] = std::min(c_ubdr - 1, my_iupper[_k]);
                        }
                    }
#ifdef _zyDEBUG
                    printf("Proc %d bdr (%d,%d,%d) ~ (%d,%d,%d)\n", my_pid,
                        bdr_ilower[0], bdr_ilower[1], bdr_ilower[2],
                        bdr_iupper[0], bdr_iupper[1], bdr_iupper[2] );
#endif
                    // 从粗往细看，自己的一个粗点对应4个细点
                    for (idx_t fngb = 0; fngb < 4; fngb ++) {
                        idx_t offset[3];
                        if (d == 0) {
                            offset[0] = __builtin_popcount(0          ); offset[1] = __builtin_popcount(fngb & 0b01); offset[2] = __builtin_popcount(fngb & 0b10);
                        } else if (d == 1) {
                            offset[0] = __builtin_popcount(fngb & 0b01); offset[1] = __builtin_popcount(0          ); offset[2] = __builtin_popcount(fngb & 0b10);
                        } else if (d == 2) {
                            offset[0] = __builtin_popcount(fngb & 0b01); offset[1] = __builtin_popcount(fngb & 0b10); offset[2] = __builtin_popcount(0          );
                        }
                        for (idx_t i2 = bdr_ilower[2]; i2 <= bdr_iupper[2]; i2++)
                        for (idx_t i1 = bdr_ilower[1]; i1 <= bdr_iupper[1]; i1++)
                        for (idx_t i0 = bdr_ilower[0]; i0 <= bdr_iupper[0]; i0++) {
                            idx_t index[3] = {i0, i1, i2}, to_index[3];
                            for (idx_t _k = 0; _k < 3; _k++) {
                                if (_k == d) to_index[_k] = ncells_per_dim - 1;// 注意此时对应加密网格的下边界
                                else         to_index[_k] = (index[_k] - (c_lbdr + 1)) * amr_stride + offset[_k];
                            }
#ifdef _zyDEBUG
                            hypre_printf("index: [%d](%d, %d, %d) - to_index: [%d](%d, %d, %d)\n",
                                        my_part, index[0], index[1], index[2],
                                        ngb_part, to_index[0], to_index[1], to_index[2]);
#endif
                            HYPRE_SStructGraphAddEntries(ssgraph, my_part, index, 0, ngb_part, to_index, 0);
                        }
                    }
                }// if upper bdr included
            }
        }
#ifdef _zyDEBUG
        MPI_Barrier(MPI_COMM_WORLD);
#else
        else
#endif
        if (my_part == 1) {// AMR加密后的网格
            idx_t ngb_part = 0;
            for (idx_t d = 2; d >= 0; d--) {// 在三个方向逐次添加连接关系
                if (my_ilower[d] == 0) {
                    assert(cart_ids[d] == 0);
                    idx_t bdr_ilower[3], bdr_iupper[3];// 确定加密网格上的边界位置
                    for (idx_t _k = 0; _k < 3; _k++) {
                        bdr_ilower[_k] = my_ilower[_k];
                        if (_k == d) bdr_iupper[_k] = my_ilower[_k];
                        else         bdr_iupper[_k] = my_iupper[_k];
                    }
#ifdef _zyDEBUG
                    printf("Proc %d bdr (%d,%d,%d) ~ (%d,%d,%d)\n", my_pid,
                        bdr_ilower[0], bdr_ilower[1], bdr_ilower[2],
                        bdr_iupper[0], bdr_iupper[1], bdr_iupper[2] );
#endif
                    // 从细往粗看，要知道自己这个细点对应哪个粗点
                    for (idx_t i2 = bdr_ilower[2]; i2 <= bdr_iupper[2]; i2++)
                    for (idx_t i1 = bdr_ilower[1]; i1 <= bdr_iupper[1]; i1++)
                    for (idx_t i0 = bdr_ilower[0]; i0 <= bdr_iupper[0]; i0++) {
                        idx_t index[3] = {i0, i1, i2}, to_index[3];
                        for (idx_t _k = 0; _k < 3; _k++) {
                            if (_k == d) to_index[_k] = c_lbdr;// 注意此时对应原网格上与加密网格的下边界
                            else         to_index[_k] = index[_k] / amr_stride + c_lbdr + 1;
                        }
#ifdef _zyDEBUG
                        hypre_printf("index: [%d](%d, %d, %d) - to_index: [%d](%d, %d, %d)\n",
                                    my_part, index[0], index[1], index[2],
                                    ngb_part, to_index[0], to_index[1], to_index[2]);
#endif
                        HYPRE_SStructGraphAddEntries(ssgraph, my_part, index, 0, ngb_part, to_index, 0);
                    }
                }
                if (my_iupper[d] == ncells_per_dim - 1) {
                    assert(cart_ids[d] == np_per_dim - 1);
                    idx_t bdr_ilower[3], bdr_iupper[3];
                    for (idx_t _k = 0; _k < 3; _k++) {
                        bdr_iupper[_k] = my_iupper[_k];
                        if (_k == d) bdr_ilower[_k] = my_iupper[_k];
                        else         bdr_ilower[_k] = my_ilower[_k];   
                    }
#ifdef _zyDEBUG
                    printf("Proc %d bdr (%d,%d,%d) ~ (%d,%d,%d)\n", my_pid,
                        bdr_ilower[0], bdr_ilower[1], bdr_ilower[2],
                        bdr_iupper[0], bdr_iupper[1], bdr_iupper[2] );
#endif
                    for (idx_t i2 = bdr_ilower[2]; i2 <= bdr_iupper[2]; i2++)
                    for (idx_t i1 = bdr_ilower[1]; i1 <= bdr_iupper[1]; i1++)
                    for (idx_t i0 = bdr_ilower[0]; i0 <= bdr_iupper[0]; i0++) {
                        idx_t index[3] = {i0, i1, i2}, to_index[3];
                        for (idx_t _k = 0; _k < 3; _k++) {
                            if (_k == d) to_index[_k] = c_ubdr;// 注意此时对应原网格上与加密网格的上边界
                            else         to_index[_k] = index[_k] / amr_stride + c_lbdr + 1;
                        }
#ifdef _zyDEBUG
                        hypre_printf("index: [%d](%d, %d, %d) - to_index: [%d](%d, %d, %d)\n",
                                    my_part, index[0], index[1], index[2],
                                    ngb_part, to_index[0], to_index[1], to_index[2]);
#endif
                        HYPRE_SStructGraphAddEntries(ssgraph, my_part, index, 0, ngb_part, to_index, 0);
                    }
                }
            }
        }

        HYPRE_SStructGraphAssemble(ssgraph);

        // 建立向量
        HYPRE_SStructVectorCreate(MPI_COMM_WORLD, ssgrid, &b);// Create an empty vector object
        HYPRE_SStructVectorCreate(MPI_COMM_WORLD, ssgrid, &x);
        HYPRE_SStructVectorCreate(MPI_COMM_WORLD, ssgrid, &y);
        HYPRE_SStructVectorCreate(MPI_COMM_WORLD, ssgrid, &x0);
        HYPRE_SStructVectorSetObjectType(b, obj_type);// Set the object type for the vectors to be the same as was already set for the matrix
        HYPRE_SStructVectorSetObjectType(x, obj_type);
        HYPRE_SStructVectorSetObjectType(y, obj_type);
        HYPRE_SStructVectorSetObjectType(x0, obj_type);
        HYPRE_SStructVectorInitialize(b);// Indicate that the vector coefficients are ready to be set
        HYPRE_SStructVectorInitialize(x);
        HYPRE_SStructVectorInitialize(y);
        HYPRE_SStructVectorInitialize(x0);
        {// 填充向量数据
            data_t* vec_buf = new data_t [local_m * local_m * local_m];
            #pragma omp parallel for collapse(2) schedule(static)
            for (idx_t i2 = my_ilower[2]; i2 <= my_iupper[2]; i2++)
            for (idx_t i1 = my_ilower[1]; i1 <= my_iupper[1]; i1++)
            for (idx_t i0 = my_ilower[0]; i0 <= my_iupper[0]; i0++) {
                const idx_t elem_id = ((i2 - my_ilower[2]) * (my_iupper[1] - my_ilower[1] + 1)
                                    +   i1 - my_ilower[1] )* (my_iupper[0] - my_ilower[0] + 1)
                                    +   i0 - my_ilower[0]  ;
                vec_buf[elem_id] = ((data_t)my_part + 1.0) * (0.0 + 1.0) * cos((data_t)(i0 + i1 + i2) / 10.0 );
                // printf("x @ Part %d (%3d,%3d,%3d) = %.10e\n", my_part, i0, i1, i2, vec_buf[elem_id]);
            }
            HYPRE_SStructVectorSetBoxValues(x, my_part, my_ilower, my_iupper, 0, vec_buf);

            // for (idx_t i = 0; i < local_m * local_m * local_m; i++) vec_buf[i] = 1.0;
            // HYPRE_SStructVectorSetBoxValues(b, my_part, my_ilower, my_iupper, 0, vec_buf);

            HYPRE_SStructVectorSetConstantValues(b, 1.0);
            HYPRE_SStructVectorSetConstantValues(y, 0.0);
            delete [] vec_buf; vec_buf = nullptr;
        }
        HYPRE_SStructVectorAssemble(b);
        HYPRE_SStructVectorAssemble(x);
        HYPRE_SStructVectorAssemble(y);// 对y也要创建
        HYPRE_SStructVectorAssemble(x0);
        HYPRE_SStructVectorCopy(x, x0);

        HYPRE_SStructMatrixCreate(MPI_COMM_WORLD, ssgraph, &A);// Create an empty matrix object
        HYPRE_SStructMatrixSetObjectType(A, obj_type);
        HYPRE_SStructMatrixInitialize(A);// Indicate that the matrix coefficients are ready to be set 
        idx_t stencil_indices[num_diag];
        for (idx_t j = 0; j < num_diag; j++)
            stencil_indices[j] = j;
        {// 填充矩阵数据
            data_t* mat_buf = new data_t [local_m * local_m * local_m * num_diag];

            #pragma omp parallel for collapse(2) schedule(static)
            for (idx_t i2 = my_ilower[2]; i2 <= my_iupper[2]; i2++)
            for (idx_t i1 = my_ilower[1]; i1 <= my_iupper[1]; i1++)
            for (idx_t i0 = my_ilower[0]; i0 <= my_iupper[0]; i0++) {
                data_t* dst_ptr = mat_buf + (((i2 - my_ilower[2]) * (my_iupper[1] - my_ilower[1] + 1)
                                          +    i1 - my_ilower[1] )* (my_iupper[0] - my_ilower[0] + 1)
                                          +    i0 - my_ilower[0]  ) * num_diag;
                dst_ptr[0] = dst_ptr[1] = dst_ptr[2] = -1.0;
                dst_ptr[3] =  6.0;
                dst_ptr[4] = dst_ptr[5] = dst_ptr[6] = -1.0;

                // Set boundary conditions
                if (i0 == 0                 ) dst_ptr[0] = 0.0;
                if (i0 == ncells_per_dim - 1) dst_ptr[6] = 0.0;
                if (i1 == 0                 ) dst_ptr[1] = 0.0;
                if (i1 == ncells_per_dim - 1) dst_ptr[5] = 0.0;
                if (i2 == 0                 ) dst_ptr[2] = 0.0;
                if (i2 == ncells_per_dim - 1) dst_ptr[4] = 0.0;

                if (my_part == 0) {// 如果本进程处于原网格
                    // Zero out dummy unknowns (原网格被加密的区域)
                    if (i0 > c_lbdr && i0 < c_ubdr && i1 > c_lbdr && i1 < c_ubdr && i2 > c_lbdr && i2 < c_ubdr) {
                        dst_ptr[3] = 1.0;
                        dst_ptr[0] = dst_ptr[1] = dst_ptr[2] = dst_ptr[4] = dst_ptr[5] = dst_ptr[6] = 0.0;
                    }
                    // Zero out off-diagonal stencil couplings at fine/coarse interface (原网格与被加密部分接壤的区域)
                    // And Fix diagonal entries at part boundaries
                    if (i1 > c_lbdr && i1 < c_ubdr && i2 > c_lbdr && i2 < c_ubdr) {
                        if (i0 == c_lbdr) { dst_ptr[6] = 0.0; dst_ptr[3] = 7.666666666666666; }
                        if (i0 == c_ubdr) { dst_ptr[0] = 0.0; dst_ptr[3] = 7.666666666666666; }
                    }
                    if (i0 > c_lbdr && i0 < c_ubdr && i2 > c_lbdr && i2 < c_ubdr) {
                        if (i1 == c_lbdr) { dst_ptr[5] = 0.0; dst_ptr[3] = 7.666666666666666; }
                        if (i1 == c_ubdr) { dst_ptr[1] = 0.0; dst_ptr[3] = 7.666666666666666; }
                    }
                    if (i0 > c_lbdr && i0 < c_ubdr && i1 > c_lbdr && i1 < c_ubdr) {
                        if (i2 == c_lbdr) { dst_ptr[4] = 0.0; dst_ptr[3] = 7.666666666666666; }
                        if (i2 == c_ubdr) { dst_ptr[2] = 0.0; dst_ptr[3] = 7.666666666666666; }
                    }
                }
                else {// 如果本进程处于细网格
                    idx_t cnt_at_fbdr = 0;
                    if (i0 ==           0       ) cnt_at_fbdr ++;
                    if (i0 == ncells_per_dim - 1) cnt_at_fbdr ++;
                    if (i1 ==           0       ) cnt_at_fbdr ++;
                    if (i1 == ncells_per_dim - 1) cnt_at_fbdr ++;
                    if (i2 ==           0       ) cnt_at_fbdr ++; 
                    if (i2 == ncells_per_dim - 1) cnt_at_fbdr ++;
                    assert(cnt_at_fbdr <= 3);
                    // if (cnt_at_fbdr > 0) printf("Proc %d (%d,%d,%d) cnt %d\n", my_pid, i0, i1, i2, cnt_at_fbdr);
                    if      (cnt_at_fbdr == 1) dst_ptr[3] = 5.666666666666666;// Fine part 2D borders (即加密网格的6个面)
                    else if (cnt_at_fbdr == 2) dst_ptr[3] = 5.333333333333333;// Fine part 1D borders (即加密网格的12条棱)
                    else if (cnt_at_fbdr == 3) dst_ptr[3] = 5.0              ;// Fine part 0D borders (即加密网格的8个顶点)
                }
            }
            HYPRE_SStructMatrixSetBoxValues(A, my_part, my_ilower, my_iupper, 0, num_diag, stencil_indices, mat_buf);
            delete [] mat_buf; mat_buf = nullptr;

            idx_t * nentries = new idx_t [local_m * local_m * local_m];
            #pragma omp parallel for schedule(static)
            for (idx_t i = 0; i < local_m * local_m * local_m; i++) nentries[i] = num_diag;
            
            // 添加上非stencil的非零元数值
            data_t offd_val = -0.666666666666666;
            if (my_part == 0) {// 原来的网格
                for (idx_t d = 2; d >= 0; d--) {// 在三个方向逐次添加非零元
                    if (my_ilower[d] <= c_lbdr && c_lbdr <= my_iupper[d]) {
                        idx_t bdr_ilower[3], bdr_iupper[3];// 确定原网格上的边界位置
                        for (idx_t _k = 0; _k < 3; _k++) {
                            if (_k == d) {
                                bdr_ilower[_k] = c_lbdr;
                                bdr_iupper[_k] = c_lbdr;
                            } else {
                                bdr_ilower[_k] = std::max(c_lbdr + 1, my_ilower[_k]);
                                bdr_iupper[_k] = std::min(c_ubdr - 1, my_iupper[_k]);
                            }
                        }
                        // 从粗往细看，自己的一个粗点对应4个细点
                        for (idx_t fngb = 0; fngb < 4; fngb ++) {
                            for (idx_t i2 = bdr_ilower[2]; i2 <= bdr_iupper[2]; i2++)
                            for (idx_t i1 = bdr_ilower[1]; i1 <= bdr_iupper[1]; i1++)
                            for (idx_t i0 = bdr_ilower[0]; i0 <= bdr_iupper[0]; i0++) {
                                idx_t index[3] = {i0, i1, i2};
                                const idx_t elem_id = ((i0 - my_ilower[0]) * local_m + i1 - my_ilower[1]) * local_m + i2 - my_ilower[2];
#ifdef _zyDEBUG
                                hypre_printf("setval: [%d](%d, %d, %d) - %d-th val %.4e\n",
                                        my_part, index[0], index[1], index[2], nentries[elem_id], offd_val);
#endif
                                HYPRE_SStructMatrixSetValues(A, my_part, index, 0, 1, &nentries[elem_id], &offd_val);
                                nentries[elem_id]++;// 计数增加
                            }
                        }
                    }// if lower bdr included
                    if (my_ilower[d] <= c_ubdr && c_ubdr <= my_iupper[d]) {
                        idx_t bdr_ilower[3], bdr_iupper[3];
                        for (idx_t _k = 0; _k < 3; _k++) {
                            if (_k == d) {
                                bdr_ilower[_k] = c_ubdr;
                                bdr_iupper[_k] = c_ubdr;
                            } else {
                                bdr_ilower[_k] = std::max(c_lbdr + 1, my_ilower[_k]);
                                bdr_iupper[_k] = std::min(c_ubdr - 1, my_iupper[_k]);
                            }
                        }
                        // 从粗往细看，自己的一个粗点对应4个细点
                        for (idx_t fngb = 0; fngb < 4; fngb ++) {
                            for (idx_t i2 = bdr_ilower[2]; i2 <= bdr_iupper[2]; i2++)
                            for (idx_t i1 = bdr_ilower[1]; i1 <= bdr_iupper[1]; i1++)
                            for (idx_t i0 = bdr_ilower[0]; i0 <= bdr_iupper[0]; i0++) {
                                idx_t index[3] = {i0, i1, i2};
                                const idx_t elem_id = ((i0 - my_ilower[0]) * local_m + i1 - my_ilower[1]) * local_m + i2 - my_ilower[2];
#ifdef _zyDEBUG
                                hypre_printf("setval: [%d](%d, %d, %d) - %d-th val %.4e\n",
                                        my_part, index[0], index[1], index[2], nentries[elem_id], offd_val);
#endif
                                HYPRE_SStructMatrixSetValues(A, my_part, index, 0, 1, &nentries[elem_id], &offd_val);
                                nentries[elem_id]++;// 计数增加
                            }
                        }
                    }// if upper bdr included
                }
            }
#ifdef _zyDEBUG
            MPI_Barrier(MPI_COMM_WORLD);
#else
            else
#endif
            if (my_part == 1) {// AMR加密后的网格
                for (idx_t d = 2; d >= 0; d--) {// 在三个方向逐次添加连接关系
                    if (my_ilower[d] == 0) {
                        assert(cart_ids[d] == 0);
                        idx_t bdr_ilower[3], bdr_iupper[3];// 确定加密网格上的边界位置
                        for (idx_t _k = 0; _k < 3; _k++) {
                            bdr_ilower[_k] = my_ilower[_k];
                            if (_k == d) bdr_iupper[_k] = my_ilower[_k];
                            else         bdr_iupper[_k] = my_iupper[_k];
                        }
                        // 从细往粗看，要知道自己这个细点对应哪个粗点
                        for (idx_t i2 = bdr_ilower[2]; i2 <= bdr_iupper[2]; i2++)
                        for (idx_t i1 = bdr_ilower[1]; i1 <= bdr_iupper[1]; i1++)
                        for (idx_t i0 = bdr_ilower[0]; i0 <= bdr_iupper[0]; i0++) {
                            idx_t index[3] = {i0, i1, i2};
                            const idx_t elem_id = ((i0 - my_ilower[0]) * local_m + i1 - my_ilower[1]) * local_m + i2 - my_ilower[2];
#ifdef _zyDEBUG
                            hypre_printf("setval: [%d](%d, %d, %d) - %d-th val %.4e\n",
                                    my_part, index[0], index[1], index[2], nentries[elem_id], offd_val);
#endif
                            HYPRE_SStructMatrixSetValues(A, my_part, index, 0, 1, &nentries[elem_id], &offd_val);
                            nentries[elem_id]++;// 计数增加
                        }
                    }
                    if (my_iupper[d] == ncells_per_dim - 1) {
                        assert(cart_ids[d] == np_per_dim - 1);
                        idx_t bdr_ilower[3], bdr_iupper[3];
                        for (idx_t _k = 0; _k < 3; _k++) {
                            bdr_iupper[_k] = my_iupper[_k];
                            if (_k == d) bdr_ilower[_k] = my_iupper[_k];
                            else         bdr_ilower[_k] = my_ilower[_k];   
                        }
                        for (idx_t i2 = bdr_ilower[2]; i2 <= bdr_iupper[2]; i2++)
                        for (idx_t i1 = bdr_ilower[1]; i1 <= bdr_iupper[1]; i1++)
                        for (idx_t i0 = bdr_ilower[0]; i0 <= bdr_iupper[0]; i0++) {
                            idx_t index[3] = {i0, i1, i2};
                            const idx_t elem_id = ((i0 - my_ilower[0]) * local_m + i1 - my_ilower[1]) * local_m + i2 - my_ilower[2];
#ifdef _zyDEBUG
                            hypre_printf("setval: [%d](%d, %d, %d) - %d-th val %.4e\n",
                                    my_part, index[0], index[1], index[2], nentries[elem_id], offd_val);
#endif
                            HYPRE_SStructMatrixSetValues(A, my_part, index, 0, 1, &nentries[elem_id], &offd_val);
                            nentries[elem_id]++;// 计数增加
                        }
                    }
                }
            }
            delete [] nentries; nentries = nullptr;
        }// fill in matrix data
        HYPRE_SStructMatrixAssemble(A);// a collective call finalizing the matrix assembly.

        for (idx_t test = 0; test < TEST_CNT; test++) {
            HYPRE_SStructVectorCopy(x0, x);
            
            data_t b_dot, x_dot, Ab_dot;
            idx_t ret;        
            // 做spmv检验一下数据传对没有
            if (obj_type == HYPRE_PARCSR) {
                HYPRE_SStructMatrixGetObject(A, (void **) &par_A);
                HYPRE_SStructVectorGetObject(b, (void **) &par_b);
                HYPRE_SStructVectorGetObject(x, (void **) &par_x);
                HYPRE_SStructVectorGetObject(y, (void **) &par_y);

                ret = HYPRE_ParVectorInnerProd(par_b, par_b, &b_dot);
                ret = HYPRE_ParVectorInnerProd(par_x, par_x, &x_dot);
                HYPRE_ParCSRMatrixMatvec(1.0, par_A, par_b, 0.0, par_y);
                ret = HYPRE_ParVectorInnerProd(par_y, par_y, &Ab_dot);

                // 用A*x作为右端项b
                HYPRE_ParCSRMatrixMatvec(1.0, par_A, par_x, 0.0, par_b);
                ret = HYPRE_ParVectorInnerProd(par_b, par_b, &b_dot);
                // 初始解设为0
                HYPRE_ParVectorSetConstantValues(par_x, 0.0);
            } else { assert(obj_type == HYPRE_SSTRUCT);
                ret = HYPRE_SStructInnerProd(b, b, &b_dot);
                ret = HYPRE_SStructInnerProd(x, x, &x_dot);
                HYPRE_SStructMatrixMatvec(1.0, A, b, 0.0, y);
                ret = HYPRE_SStructInnerProd(y, y, &Ab_dot);

                // 用A*x作为右端项b
                HYPRE_SStructMatrixMatvec(1.0, A, x, 0.0, b);
                ret = HYPRE_SStructInnerProd(b, b, &b_dot);
                // 初始解设为0
                HYPRE_SStructVectorSetConstantValues(x, 0.0);
            }
            
            if (my_pid == 0) {
                printf("(  b,   b) = %.20e\n",  (double)b_dot);
                printf("(  x,   x) = %.20e\n",  (double)x_dot);
                printf("(A*b, A*b) = %.20e\n", (double)Ab_dot);
            }

            buildup_solver(solver_name, prec_name, config);

            idx_t num_iterations = 0;
            data_t t_setup = 0.0, t_solve = 0.0, final_res_norm = 0.0;
            
            if (obj_type == HYPRE_SSTRUCT)
                setup_and_solve(solver_name, prec_name, A, b, x, t_setup, t_solve, final_res_norm, num_iterations);
            else if (obj_type == HYPRE_PARCSR) {
                setup_and_solve(solver_name, prec_name, par_A, par_b, par_x, t_setup, t_solve, final_res_norm, num_iterations);
                HYPRE_SStructVectorGather(b);
                HYPRE_SStructVectorGather(x);
            }

            // 计算真实残差
            data_t r_nrm2 = 0.0, b_nrm2 = 0.0;
            check_residual(prec_name, A, x, b, y, r_nrm2, b_nrm2);
            if (my_pid == 0) {
                printf("\033[1;35mtrue ||r|| = %20.16e ||r||/||b||= %20.16e\033[0m\n", r_nrm2, r_nrm2/b_nrm2);
                printf("Iterations = %d\n", num_iterations);
                printf("Time cost %.5f %.5f %.5f %d\n", (double)t_setup, (double)t_solve, (double)(t_setup + t_solve), num_iterations);
                printf("Final Relative Residual Norm = %e\n", (double)final_res_norm);
                printf("\n");
            }

            records[test].iter = num_iterations;
            records[test].setup = t_setup;
            records[test].solve = t_solve;

            destroy_solver(solver_name, prec_name);
        }// test loop

        if (my_pid == 0) {
            std::sort(records, records + TEST_CNT);
            const TEST_RECORD & best = records[0], & worst = records[TEST_CNT - 1];
            TEST_RECORD avg;
            for (idx_t i = 0; i < TEST_CNT; i++) {
                avg.setup += records[i].setup;
                avg.solve += records[i].solve;
                avg.prec  += records[i].prec;
                avg.iter  += records[i].iter;
            }
            avg.setup /= TEST_CNT;
            avg.solve /= TEST_CNT;
            avg.prec  /= TEST_CNT;
            avg.iter  /= TEST_CNT;
            data_t min_tot = best .setup + best .solve;
            data_t max_tot = worst.setup + worst.solve;
            data_t avg_tot = avg  .setup + avg  .solve;
            printf("\n  Summary of %d tests:\n", TEST_CNT);
            printf("     Setup time  Solve time  #Iter  Total time\n");
            printf("Min    %.6f    %.6f  %5d  %.6f\n", best .setup, best .solve, best .iter, min_tot);
            printf("Avg    %.6f    %.6f  %5d  %.6f\n", avg  .setup, avg  .solve, avg  .iter, avg_tot);
            printf("Max    %.6f    %.6f  %5d  %.6f\n", worst.setup, worst.solve, worst.iter, max_tot);    
        }
    
        HYPRE_SStructMatrixDestroy(A);
        HYPRE_SStructVectorDestroy(b); HYPRE_SStructVectorDestroy(x); HYPRE_SStructVectorDestroy(y); HYPRE_SStructVectorDestroy(x0); 
        HYPRE_SStructGraphDestroy(ssgraph);
        HYPRE_SStructStencilDestroy(stencil);
        HYPRE_SStructGridDestroy(ssgrid);
    }
    HYPRE_Finalize();
    MPI_Finalize();
    return 0;
}